using System;
using System.Reflection;
using NBody;
using NBody.InitialConditions;
using OptionParser;
using Gsl;

[assembly: AssemblyTitle ("NBodyFromDF")]
[assembly: AssemblyVersion ("1.0.*")]
[assembly: AssemblyDescription ("Build a N-body system with a specified DF")]
[assembly: AssemblyCopyright ("2005 Joseph D. MacMillan")]

public class FromDFOptions : OutputGetOptions
{
    [Option ("-n", "Number of particles")]
    public int N;
    
    [Option ("-c", "Adjust for centre of mass")]
    public bool CentreOfMass;
    
    public FromDFOptions()
    {
        N = 10000;
    }
}

// can't use Sphere.GenerateFromDF because of radial-only orbits
public class NBodyFromDF
{
    double total_mass;   
    double beta;
    double a;
    
    public NBodyFromDF(double beta, double a)
    {
        this.beta = beta;
        this.a = a;
        total_mass = 4.0 * Math.PI * beta * a;
    }
    
    public static void Main(string[] args)
    {
        FromDFOptions opts = new FromDFOptions();
        opts.ProcessArgs(args);
        
        NBodySystem s = new NBodyFromDF(0.25/Math.PI, 1.0).GenerateFromDF(opts.N, 0.0, 2);
        
        if (opts.CentreOfMass)
            s.AdjustForCentreOfMass();
        
        s.Write(opts.OutFile);
    }
    
    private NBodySystem GenerateFromDF(int n, double time, int seed)
    {
        NBodySystem sys = new NBodySystem(n, time);
        Random urandom;
        Random grandom;
        if (seed > 0)
        {
            urandom = new Random(seed);
            grandom = new RandomGaussian(seed);
        }
        else
        {
            urandom = new Random();
            grandom = new RandomGaussian();
        }
        
        double E0 = -4.0 * Math.PI * beta;
        double mass = total_mass / (double)n;
        
        for (int i = 0; i < n; i++)
        {
            double radius, E, w1;
            do
            {
                double s = Math.Sqrt(0.5) * grandom.NextDouble(); // check
                E = E0 * (1.0 + s * s);
                w1 = 2.0 * Math.PI * urandom.NextDouble();
                double w1_plus;
                if (w1 > Math.PI)
                    w1_plus = 2.0 * Math.PI - w1;
                else
                    w1_plus = w1;
            
                double v = InverseErf(1.0 - w1_plus / Math.PI);
                radius = a * Math.Exp((v * v - E + E0) / E0);
            }
            while (radius == 0.0);
            
            double x, y, z;
            PickShell(urandom, radius, out x, out y, out z);
                        
            double velocity = Math.Sqrt(2.0 * E + 2.0 * E0 * Math.Log(radius / a) - 2.0 * E0);
                        
            if (w1 > Math.PI)
                velocity = -velocity;
            
            double vx = velocity * x / radius;
            double vy = velocity * y / radius;
            double vz = velocity * z / radius;
            
            sys.AddParticle(new Particle((double)mass, (double)x, (double)y, (double)z, (double)vx, (double)vy, (double)vz));
        }
        
        return sys;
    }
    
    // uses binary searching to find inverse
    const double INF = 10000.0;
    const double tolerance = 1e-7;
    private static double InverseErf(double x)
    {
        if (x == 0.0) return 0.0;
        else if (x < -1.0 || x > 1.0)
            throw new ApplicationException("InverseErf: x = " + x);
         
        double min = -INF;
        double max = INF;
                 
        double error = INF;
        double midpoint = 0.0;
        while (error > tolerance)
        {
            midpoint = min + 0.5 * (max - min);
            double f = SpecialFunctions.Erf(midpoint);
            if (f > x )
                max = midpoint;
            else
                min = midpoint;
            error = Math.Abs((f - x) / x);
        }
         
        return midpoint;
    }
    
    private NBodySystem GenerateFromDF2(int n, double time, int seed)
    {
        NBodySystem s = new NBodySystem(n, time);
        Random random;
        if (seed > 0)
            random = new Random(seed);
        else
            random = new Random();
        
        double m_frac = 1.0 * total_mass;
        double g_max = 100.0; // this depends on g_function
                                
        for (int i = 0; i < n; i++)
        {
            double mass = total_mass / (double)n;
            double radius = mass_inverse(m_frac * random.NextDouble());
            
            double x, y, z;
            PickShell(random, radius, out x, out y, out z);
            
            double q, test;
            do
            {
                q = random.NextDouble();
                test = g_max * random.NextDouble();
            } while (test > g(q));
            
            double velocity = q * Math.Sqrt(-2.0 * phi(radius));
            
            // velocity must be radial, either inward or out
            double direction = random.NextDouble();
            if (direction < 0.5) velocity = -1.0 * velocity;
            
            double vx = velocity * x / radius;
            double vy = velocity * y / radius;
            double vz = velocity * z / radius;
            
            s.AddParticle(new Particle((double)mass, (double)x, (double)y, (double)z, (double)vx, (double)vy, (double)vz));
        }
        
        return s;
    }
    
    private double mass_inverse(double m)
    {
        return m / 4.0 / Math.PI / beta;
    }
    
    private double phi(double r)
    {
        return 4.0 * Math.PI * Math.Log(r / a) - 4.0 * Math.PI;
    }
    
    private double g(double v)
    {
        return v * v / Math.Sqrt(1.0 - v * v);
    }
    
    private void PickShell(Random random, double radius, out double x, out double y, out double z)
    {
        double theta = Math.Acos(1.0 - 2.0 * random.NextDouble());
        double phi = 2.0 * Math.PI * random.NextDouble();
    
        x = radius * Math.Cos(phi) * Math.Sin(theta);
        y = radius * Math.Sin(phi) * Math.Sin(theta);
        z = radius * Math.Cos(theta);
    }
}
